import numpy as np
import pandas as pd
from numpy import diff
from scipy.signal import savgol_filter
from scipy.integrate import simps
from tkinter import messagebox

#############################################################################################################
# General API

def find_Nearest_X_Point(x,y, point_Y_value ):   #传入的xy是一个数组
    nearest_X_point=-50
    for i in range(len(y)-1):
            value1_pre=y[i]-point_Y_value
            value2_pre=y[i+1]-point_Y_value
            if value1_pre*value2_pre<0:
               if abs(value1_pre)<abs(value2_pre):
                  nearest_X_point=x[i]
                  break
               else:
                  nearest_X_point=x[i+1]
                  break
    return nearest_X_point

def Find_Nearest_Value(data , value):   #找到data数组中距离value最近的值的索引
    array=np.array(data)
    idx = (np.abs(array-value)).argmin()
    #value_t = round(float(array[idx]),2)
    return idx

#############################################################################################################
#############################################################################################################
# 3-Port and 2-Port IV Extraction Algorithm

def Smooth_gm_Transconductance(gm,width):
    for i in range(len(gm)-1):
       if gm[i]>=0.7*(gm[i-1]+gm[i+1]) and i>=1 and gm[i-1]/(width*0.001)>2E-3:    ##认为2mS/mm是gm已经跨过了快速增长的代码陷阱区。
          if gm[i+1]<=2.5*gm[i]:
             if gm[i+1]>=0.7*gm[i]:
                gm[i]=0.5*(gm[i-1]+gm[i+1])
             else:
                gm[i+1]=0.85*gm[i]
                gm[i]=0.5*(gm[i-1]+gm[i+1])
          else:
             gm[i+1]=1.5*gm[i]
             gm[i]=0.5*(gm[i-1]+gm[i+1])
       if abs(gm[i])<=0.375*(gm[i-1]+gm[i+1]) and i>=1 and gm[i-1]/(width*0.001)>2E-3:
          if gm[i+1]<=2.5*gm[i]:
             if gm[i+1]>=0.7*gm[i]:
                gm[i]=0.5*(gm[i-1]+gm[i+1])
             else:
                gm[i+1]=0.85*gm[i]
                gm[i]=0.5*(gm[i-1]+gm[i+1])
          else:
             gm[i+1]=1.5*gm[i]
             gm[i]=0.5*(gm[i-1]+gm[i+1])
#gm的Savgol滤波（自动平滑）算法：
    gm_smooth = savgol_filter(gm, 21, 2, mode='interp')  #【重要】Savgol滤波法 平滑处理gm图像，参数21是窗口长度，参数2是多项式阶数.根据图像处理效果调节
  #注意，smooth算法的nearest与插值算法中的nearest不同。对于后者，这种算法非常粗糙，容易引入锯齿以及图形形状的失真。linear或者spline会好很多。
  #但smooth算法是一种卷积运算，采用nearest做mode效果还可以。不过综合对比之后，还是interp mode得出的GVS会偏大一些，图像的锯齿也更少。
    return gm_smooth

def Extract_Vth_igleak_vgleak(vg,ig,id,width):

   Vth:float=-50
   ig_leak:float=0
   vg_leak_max:float=-50

#提取Vth及负栅压时最大栅漏电ig_leak
   for i in range(len(id)):    #以下的所有i-1都可以通过在for循环声明里table.nrows-1来一劳永逸。因为表头这一行不能计入数据采集array。
        if(id[i]>1E-5*width*0.001):   #单位是A，1E-6对应1uA，则1E-5对应10uA；进行归一化得10uA/mm，认定此时器件开启，记录对应的Vth.
           Vth=vg[i]
           break
   for i in range(len(ig)):
     if isinstance(ig[i],float):
        if vg[i]<0 and ig_leak<abs(ig[i]):  #负栅压条件判断 + 手动找到最大栅漏电的算法
            ig_leak=abs(ig[i])
            vg_leak_max=vg[i]
   return Vth,ig_leak,vg_leak_max

def Extract_gmax_Vgmax(gm_smooth,vg):
   g_max:float=0
   vg_max:float=-50
#提取gmax
#寻找特定Vd下，平滑转移曲线上最大的跨导gmax（gm_smooth已经过平滑处理，排除高频噪声尖峰的干扰），定义为gmax.同时，抓取gmax对应的Vg_max点.
   for i in range(len(gm_smooth)):
        if(g_max<gm_smooth[i]):   #i-2是为了平衡gm_smooth的第一项为excel表中的第三项，避免抓取失败
           g_max=gm_smooth[i]
           vg_max=vg[i+1]
   return g_max,vg_max

def Extract_GVS_Vg_Boundary(gm_smooth,g_max,vg_max,vg):
#提取GVS
   gmk=0.8*g_max        #GVS的临界定义：80%的gmax对应的边界。
   vk1:float=-50  #设置为-50是为了检验程序运行是否出错.
   vk2:float=-50  #设置为-50是为了检验程序运行是否出错.
   get_GVS:int=0
   for i in range(len(gm_smooth)):
      if i>=1:
          gm1_pre=gm_smooth[i-1]-gmk
          gm2_pre=gm_smooth[i]-gmk
          if gm1_pre*gm2_pre<0 and get_GVS!=1:
             if abs(gm1_pre)<abs(gm2_pre):
                vk_pre=vg[i-1]  #table.cell_value(i-2,col_vg)
             else:
                vk_pre=vg[i]  #table.cell_value(i-1,col_vg)
             if vk_pre<vg_max: #and abs(gm1_pre)<=0.05*gmk and abs(gm2_pre)<=0.05*gmk:
                  vk1=vk_pre
             else:
                  vk2=vk_pre
                  get_GVS=1
                  break
   if(get_GVS==0):
     vk2=vg[len(gm_smooth)]    #table.cell_value(len(gm_smooth)+1,col_vg)
     #如果在整个sweep周期内vk2都没有被赋值，说明在测量范围内未找到gm右侧降到80%GVS的点。这时将右端赋值为测量VG的【右边缘】，一般也就是测试的上限。
   GVS=vk2-vk1
   return GVS,vk1,vk2


def Extract_Idss_Ron_Vk(vd,id,width):  #Vk提取，如果对应的Idss归一化后在100mA/mm以下，则降低i_k标准，防止从未开启时的输出曲线抓取Vk
    i_dss=0
    for i in range(len(id)):
        if isinstance (id[i],float):
          if(i_dss<id[i]):
             i_dss=id[i]   #寻找特定栅压条件下，输出曲线上最大的电流，定义为Idss.
    Ron:float=0
    vd_id_zero=find_Nearest_X_Point(vd,id,0)
    if vd_id_zero==-50: vd_id_zero=0
    for i in range(len(id)):
       if vd[i]>=vd_id_zero and vd[i]>=0.03:  #后面一个判据是防止用栅漏电部分(Id<0)数据计算. 阈值定在30mV避免pulse系统的噪声干扰。
         if Ron==0 or ((vd[i]-vd_id_zero)/id[i])<=1.001*Ron:  #寻找输出曲线过原点割线斜率最大、即最小的Ron
            Ron=(vd[i]-vd_id_zero)/id[i]
            # print("Ron,vd[i],vd_id_zero,id[i]",Ron,vd[i],vd_id_zero,id[i])
    if i_dss*1000*1000/width>100:
       Vk=find_Nearest_X_Point(vd,id,0.9*i_dss)
    elif i_dss*1000*1000/width>50:
       Vk=find_Nearest_X_Point(vd,id,0.7*i_dss)
    else:
       Vk=-10
    return i_dss,Ron,Vk


def Extract_VBD_or_Vmea(vd,id,id_bd_th):    #id_bd_th是反归一化之后的结果
    vd_breakdown:float=-50
    vmax_mea:float=10
    BD:int=0
    for i in range(len(id)):
        if id[i]>id_bd_th and BD==0 :   #击穿发生与否，判断id是否达到击穿阈值id_bd_th(在控制台指定)
              vd_breakdown=vd[i]
              BD=1  #记录此时击穿已发生。防止重复赋值vd_breakdown
        if vd[i]>vmax_mea:
              vmax_mea=vd[i]  #定出测试范围：测试过程中出现的最大漏压
    if BD==0: vd_breakdown=vmax_mea
    return vd_breakdown,vmax_mea

#############################################################################################################
#############################################################################################################
# CV Extraction Algorithm

def Smooth_Capacitance_data(filename, vb , cap_raw , **kwargs):   #滤波，剔除，再滤波
    index_start=0
    index_stop=len(cap_raw)
    if vb[0]<vb[1]:
        mode="mi"   #monotonically increasing 默认单调递增，也可以输入'md'，单调递减
    else:
        mode="md"
    for k,v in kwargs.items():
        if k=="index_start":
            index_start=v
        if k=="index_stop":
            index_stop=v
        if k=="mode":
            mode=v

    if isinstance(vb,np.ndarray):   vb=vb.tolist()
    if isinstance(cap_raw,np.ndarray):   cap_raw=cap_raw.tolist()
    vb0 = vb[0:index_start]
    cap0 = cap_raw[0:index_start]
    vb1 = vb[index_start:index_stop]
    vb2 = vb[index_stop:len(vb)]
    cap1 = cap_raw[index_start:index_stop]
    cap2 = cap_raw[index_stop:len(cap_raw)]

    cap1 = savgol_filter(cap1 , 5, 2, mode='interp')
    vb1_trimmed=[vb1[0]]
    cap1_trimmed=[cap1[0]]
    vb1_deleted=[]
    cap1_deleted=[]
    for i in range(1,len(cap1),1):
        if mode=="mi":
            if cap1[i]>=cap1_trimmed[len(cap1_trimmed)-1] and vb1[i]>vb1_trimmed[len(vb1_trimmed)-1]:   #只往cap_trimmed数组里放数据。9.20变严格单增(>)为单增(>=)
                vb1_trimmed.append(vb1[i])
                cap1_trimmed.append(cap1[i])
            else:
                vb1_deleted.append(vb1[i])
                cap1_deleted.append(cap1[i])
        elif mode=="md":
            if cap1[i]<=cap1_trimmed[len(cap1_trimmed)-1] and vb1[i]>vb1_trimmed[len(vb1_trimmed)-1]:   #只往cap_trimmed数组里放数据。9.20变严格单减(<)为单减(<=)
                vb1_trimmed.append(vb1[i])
                cap1_trimmed.append(cap1[i])
            else:
                vb1_deleted.append(vb1[i])
                cap1_deleted.append(cap1[i])
    try:
        cap1_trimmed=savgol_filter(cap1_trimmed, 5, 2, mode='interp')
    except:
        # print("cap1",cap1)
        df_cap=pd.DataFrame(cap1)
        df_cap.to_clipboard(excel=True, index=False)
        # print("vb1_trimmed",vb1_trimmed)
        # print("cap1_trimmed",cap1_trimmed)
        vb1_trimmed=vb1
        cap1_trimmed=cap1
        inform_warning="您指定的数据文件：",filename,"不适用于强力滤波。","已自动调整为常规滤波算法。若要继续采用强力滤波算法，请剪裁数据文件后重新导入。"
        inform_warning='\n'.join(inform_warning)
        messagebox.showwarning('提示',inform_warning)

    cap_sum=cap0
    cap_sum.extend(cap1_trimmed)
    cap_sum.extend(cap2)
    vb_sum=vb0
    vb_sum.extend(vb1_trimmed)
    vb_sum.extend(vb2)
    return vb_sum,cap_sum

def Sort_nx_d ( d , nx):            #按照d从小到大的顺序重新排列d和nx
    d0 = sorted(d)
    d_new=[]
    nx_new=[]
    for i in range(len(d0)):
        index_sort = d.index(d0[i])
        d_new.append(d0[i])
        nx_new.append(nx[index_sort])
    return d_new , nx_new

def Extract_Vth_from_CV (filename, vb_raw , c_raw , mode_filter="Normal") :
      if mode_filter == "Normal":
         vb = vb_raw
         c = savgol_filter(c_raw,5,2,mode='interp')
      if mode_filter == "Raw":
         vb = vb_raw
         c = c_raw
      if mode_filter == "Smooth":
         vb , c = Smooth_Capacitance_data (filename, vb_raw , c_raw)
      c_max = max(c)
      vth = find_Nearest_X_Point(vb , c , 0.5*c_max)
      vth = round(vth,2)
      return vth

def Extract_Nsh_from_CV_Method1(vb , cap_raw ,A):   #基于CdV/Aq直接积分法
    cap_filter=True

    q: float = 1.6E-19  # 元电荷，单位C.
    if cap_filter:
        cap = savgol_filter(cap_raw, 5, 2, mode='nearest')
    else:
        cap = cap_raw
    Qsh_1: float = (simps(cap, vb)) / A
    nsh_1=Qsh_1/q
    nsh_1=abs(nsh_1*1E-4)  #将m-2转化成cm-2单位，并将负号调正.
    nsh_1_printed='%.2e' % nsh_1
    #print('按照直接积分算法，该CV环结构对应的2DEG面密度为: ',nsh_1_printed ,' cm-2')
    return nsh_1_printed

def Extract_Nsh_d_nx_from_CV_Method2(filename, vb , cap_raw , A ,k , mode="Normal"):  #nx-d间接积分法
    # 采用C^3算法计算2DEG面密度，同步导出体浓度nx - 厚度d的曲线关系图。
    # 注意，采用这种算法，要获得较为精确的nx-d关系，必须在CV曲线的上升段 减小偏压步进的step，提高2DEG开启过程的测量分辨率，不然d的取样离散化程度很大，各点之间的平滑过渡很难保证.
    # 另外，由于d近端的nx会在很大程度上受平滑处理算法、dV/dC的尖峰影响，造成起始点的nx偏高，从而离散数组积分得到的nsh会偏大。
    q: float = 1.6E-19  # 元电荷，单位C.
    d=[]   #d的单位是m.
    nx=[]   #nx的单位是/m^3.
    if mode=="Normal":
        cap = savgol_filter(cap_raw, 5, 2, mode='nearest')
    elif mode=="Smooth":
        vb , cap = Smooth_Capacitance_data(filename, vb,cap_raw)
    elif mode=="Raw":
        cap = cap_raw
    dvdc=diff(vb)/diff(cap)

    for i in range(len(vb)-1):
           d.append((k*A/(abs(cap[i])))  ) #d[i]=k*A/(c[i])  【】abs是否使用有待商榷.  d是空数组，采用insert注入数据.
           nx.append(abs((pow(cap[i],3)*(dvdc[i])/(pow(A,2)*k*q))) )  #【】abs是否使用有待商榷.

    d_eff=[]
    nx_eff=[]
    for i in range(len(vb)-1):
        if d[i] < 200E-9:  # 积分计算只计算200nm厚度以内的nx-d，因为对于常规异质结结构，2DEG最多分布在该范围内。不然积分值会偏大
            d_eff.append(d[i])
            nx_eff.append(nx[i])
    if len(d_eff) != 0:
        nsh=simps(nx_eff,d_eff)    #对离散数组积分求出2DEG面密度，这里取单位转化前的值进行计算，防止单位出错.
    else:
        nsh=simps(nx,d)

    for i in range(len(d)):
       d[i]=d[i]*1E9  #将d的m单位变化到nm单位
       nx[i]=nx[i]*1E-6  #将nx的m-3单位变化到cm-3单位
       #if(d[i]<=1000):   #只显示1000nm以内的数据，因为只有这部分靠近表层，成为导电层贡献。更深处可能是buffer层导电.
          #print(d[i],nx[i])

    #nx = savgol_filter(nx, 51, 2, mode= 'nearest')
    #显示积分处理的面密度结果（显示在作图之前，防止吞数据）
    nsh=abs(nsh*1E-4)  #将m-2转化成cm-2单位，并将负号调正.
    nsh_printed : str = '%.2e' % nsh
    #print('按照nx-d算法，该CV环结构对应的2DEG面密度为: ', nsh_printed ,' cm-2')
    return nsh_printed, d, nx


